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We numerically simulate turbulent Taylor-Couette flow for independently rotating in- 
ner and outer cylinders, focusing on the analogy with turbulent Rayleigh-Bcnard flow. 
Reynolds numbers of Re t — 8 • 10 3 and Re D — ±4 • 10 3 of the inner and outer cylinders, 
respectively, are reached, corresponding to Taylor numbers Ta up to 10 8 . Effective scal- 
ing laws for the torque and other system responses are found. Recent experiments with 
the Twcntc turbulent Taylor- Couette (T 3 C) setup at very high Reynolds numbers have 
revealed an optimum transport at a certain non-zero rotation rate ratio a = —uj /uji 
that depends on Ta. For large enough Ta in the numerically accessible range we find 
such an optimum at non-zero counter- rotation also in the numerics. We furthermore nu- 
merically calculate the corresponding angular velocity profiles and visualize the different 
flow structures for the various regimes. By writing the equations in a frame co-rotating 
with the outer cylinder a link is found between the local angular velocity profiles and the 
global transport quantities. 
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1. Introduction 

Taylor-Couette (TC) flow, i.e. the flow in the gap between two independently rotating 
coaxial cylinders, is among the most investigated problems in fluid mec hanics, due to its 



Or 



concep tional simplicity and to applications in process technology, see e.g. lHaim fe Pismen 



(1994). Traditionally, the driving of this system is expressed by the Reynolds numbers 



of the inner and outer cylinders, defined by Rei — r^id/v and Re Q — r ui d/v, where 
d = r ~ Ti is the gap width and v the kinematic viscosity. The geometry of a TC system 
is expressed by the radius ratio T) = ri/r and the aspect ratio V = L/d, see figure [TJ 
It was shown by Eckhardt, Grossmann & Lohse 2007 (from now on referred to as EGL 
2007) that TC flow has many similarities to Rayleigh-Benard (RB) convection, i.e. the 
thermal flow in a fluid layer heated from below and cooled from above, which will be 
discussed in detail below. 

Both RB and TC flows have been popular playgrounds for the developm ent of new con- 



cepts in fluid dynamics. Both systems have been used to study instabilities (iPfister fc Rehbere 
198lUPfister et a/.lll988llChandrasekharlll98ll:lDrazin fc ReiJl98lUBusselll967l). nonlin - 



ear dynamics and chaos (jLorenzl Il963t lAhlerd Il974t iBehrrngerl Il985t IStrogatzl Il99l 
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Figure 1: Geometry of the Rayleigh-Benard (left) and the Taylor-Couette systems (right). 
The RB system consists of two cylindrical plates, a hot one at the bottom and a cold one 
at the top of diameter D separated by a distance L. The top plate is at a temperature 
To and the bottom plate is at a temperature To + A, with A the temperature difference 
between the plates. The TC system consists of two coaxial cylinders of length L. The 
inner cylinder has the radius Ti and the angular velocity Ui , while the outer cylinder has 
the radius r Q and the angular velocity uj . 



pattern formation (Cross & Hohenbers 


1993: Bodenschatz et all 20001). and turbulence 


(ISitredall994l Grossmann & Lohsd200(l 


KadanoflteOOll LathroD et aZ.ll992&: Ahlers et al. 


2009: Lohse & Xia!l2010f). The reasons that RB and TC are so popular include: (i) These 



systems are mathematically well-defined by the Navier-Stokes equations and the appro- 
priate boundary conditions; (ii) these are closed system and thus exact global balance 
relations between the driving and the dissipation can be derived; and (iii) they are exper- 
imentally and numerically accessible with high precision, thanks to the simple geometries 
and high symmetries. 

The analogy between TC and RB may be better seen from the exact relations (EGL 
2007) between the transport quantities and the energy dissipation rates. For RB flow the 
conserved quantity that is transported is the thermal flux J = (u z 9) A.t — nd z (0) A,t of the 
temperature field where the first term is the convective contribution (u z is the vertical 
fluid velocity component) and the second term is the diffusive contribution. Here (...) A,t 
indicates the averaging over time and a horizontal plane. In the state with lowest thermal 
driving there is not yet convection. Therefore J = Jq — kAL^ 1 and the corresponding 
dissipation rate is e Ui o = since u = 0. - In TC flow, the conserved transport quantity, 
which is transported from the inner to the outer cylinder (or vice versa) is the flux 
ju) _ r 3 ((u r w)A,t — vd r (to)A,t) of the angular velocity field u, where the first term is the 
convective contribution with u r as the radial fluid velocity component and the second 
term is the diffusive contribution, cf. EGL 2007. In this case (...)a,* indicates averaging 
over time and a cylindrical surface with constant radial distance r from the axis. In the 
state with lowest driving induced by the rotating cylinders and neglecting plate effects 
from the upper and lower plates (achieved in the DNS by periodic boundary conditions 
in axial direction), the flow is laminar and purely azimuthal, ug(r) = Ar + B/r, while 
u r = u z = 0. This flow provides an angular velocity current Jtf (called Jf am in EGL2007) 
and a nonzero dissipation rate, see eqs. (|1.8|) and (|1.16|) in table 1. 

The analogy between RB and TC (EGL 2007) is highlighted when the driving in TC is 
expressed in terms of the Taylor number Ta and the angular velocity ratio a — —u> /uji of 
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Figure 2: Control parameter phase space which was numerically explored in this paper: 
(a) classical representation (Re ,Rei) and (b) (Ta,a) representation with a = —uj /uji 
and (c) (Ta, l/Ro) representation with Ro defined in eq. (I2.3[) . We fixed r\ — 5/7, T = 2tt 
and employed axial periodicity. The same color code, denoting the Taylor number, is 
maintained throughout the paper. The grey-shaded area is the Rayleigh-stable regime. 



the cylinders, while the response is given by the dimensionless transport current density 
J u divided by the corresponding molecular current density of the angular velocity from 
the inner to the outer cylinder, called the "w-Nusselt number" Nu^,. The Taylor number 
Ta is defined as Ta = a(r a — n) 2 (r + ri) 2 (uj — uii) 2 /(Av 2 ), or 



Here 



Ta=(!*d*/ry i v*){u -u i y 



4 / 4 



(1.1) 



(1.2) 

with r a = (r a + r^)/2 the arithmetic and r g = yjr^rl the geometric mean radii. are 
the angular velocities of the outer and inner cylinders, respectively; see also table [T] for 

definitions and relations. 

TC flow has been extensively investigated experimentally (IWendt 19331: Taylor 193€ ; 
Smith fc Townsendlll982l: lAndereck aLlll986|:lTong et aLlll990HLathro'p" et al. 1992&. a; 
Lewis fe Swinne vl ll999l:lvan Cils et aLll2011fel U: lHuisman et a,l\\201$\ at low and high Ta 
numbers for different ratios of the rotation frequencies a = —uj /uji 7 see the phase dia- 
gram in figure[3] However, up to now most numeri cal simulations of TC flow have been re- 
stricted to the case of pure inner cylinder rotat ion ( Fasel fc Boozll984tlCoughlin fc Marcu; 
1996tlDonal2007l . r2008l : |Pirro fc Quadriol2008j) . or e igenvalue study dGebhardt fc Grossman 
1993), o r counter- rotation at fixed a (Dong 2008). Rec ent experiments (jvan Gils et al 



inn 



2011al f6: Paoletti fc Lathrop 2011 : Huisman et al. 2012) have shown that at fixed Ta an 



optimal transport is obtained at non-zero a. 

In this paper we use direct numerical simulations (DNS) to study the influence of the 
choice of a on the flow structures and the corresponding transported angular velocity 
flux in detail for Ta numbers up to Ta = 10 8 . Figure [5] shows the cases which have been 
simulated here in the (Ta,a) as well as in the (Re , Re{) parameter spaces. Note that a 
higher density of points has been used in places where the response (Nu u , Re w ) shows 
more variation. All points have been simulated for fixed T = 2n and n — 5/7 since these 
are very similar to the parameters of the T 3 C setup. There is a significant difference, 
however, as numerically we take periodic boundary conditions in the axial direction, 
while the T 3 C system is a closed system with solid boundaries at top and bottom. 

In chapter 2 we start with a description of the numerical method that has been used. 
In chapter 3 we will discuss the validation and resolution tests that have been performed. 
In chapter 4 the global response, in terms of Nu u and the wind Reynolds number Re w , 
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Table 1: Analogous relations between RB and TC flow, leading to th e same effective scaling 
laws as derived by Eckhardt, Grossmann & Lohse (Eckhardt et al. 2007). In RB flow, the dimen- 
sionless control parameters are the Rayleigh number Ra = f3gAL 3 /(vk), the Prandtl number 
Pr = v/k, and the aspect ratio V — D/L, where A is the temperature difference between the 
cold top and hot bottom, f3 the thermal expansion coefficient, g the gravitational acceleration, 
and k the thermal diffusivity, see figure [1] The response of the system is the heat flux from 
the bottom to the top in terms of the molecular heat flux, known as the Nusselt number Nu. 
In analogy, for TC flow we define a Nusselt number Nu^ as ratio of the total and the purely 
azimuthal and laminar angular velocity flow. e Uj o is the dissipation in the purely diffusive state, 
equal to zero in RB flow, since the fluid velocity is zero and there is molecular heat transport 
only, while in TC flow e Ui o is the purely azimuthal and laminar flow dissipation rate. 



Rayleigh-Benard 



Taylor-Couette 



Conserved: temperature flux 

J = {u z 0) A ,t - Kd z {6)A,t 
Dimensionless transport: 



Nu - 



J 



Jo 

Jo = acAL -1 



Driven by: 



Ra = 



(3gAL 3 



Exact relation: 



e u — £ u — e-u.o 



Prandtl number 
Scaling: 



Pr — v Ik 
Nu oc Ra? 



(1.3) 

(1.5) 
(1.7) 

(1.9) 
(1.11) 



= (Nu — 1) Ra Pr' 2 (1.13) 
eu,o = (1.15) 



(1.17) 
(1.19) 



Conserved: angular velocity flux 

J" — r3 ((<*rW)i,( - vd r (u))A,t) (1-4) 

Dimensionless transport: 



Nu u 



2nLpJ^ 



(1.6) 



2rf rl toi - u 
J o = V Z — ~ j ( L8 ) 



Ti + r a d 
Driven by: 

1 a(r -ri) 2 (ri+r ) 2 (iJi-u} ) 2 



Ta - 

Exact relation: 
~/ ~ 



(1.10) 
(1.12) 



(Nuu) — 1) Ta a (1.14) 



d r t r a (uii - uj 



V" 



d 



Pseudo 'Prandtl' number: 

<x=(l + ^7f4: n;XJ 



Scaling: 



' o 



• o 



Nuu, oc Ta 7 



(1.16) 

(1.18) 
(1.20) 



as functions of the angular velocity ratio a will be discussed. In order to understand 
the global system response we will analyze the coherent structures in chapter 5 and 
the boundary layer profiles in chapter 6, i.e. quantities that are difficult to analyze in 
experiments. This allows us to rationalize the position of the maximum in Nu LJ (a). We 
c onclude with a brief discus s ion an d outlook to future work in chapter 7. 

Brauckmann fc Eckhardt ( 20121 ) offer a complementary direct numerical simulation 



of turbulent Taylor-Couette flow: They employ a spectral code with periodic boundary 
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Figure 3: Explored phase space (Re a , Ra) of TC flow with indepe n dently rotat 
ing inner and outer cyl i nders. Both 



experimental data (IWendtl Il933t iTavlorl Il936| 
3Hll986t iTong et ami990t lLathrop et a/Jll992al 



Smith fc Townsend] Il982t [Andereck et aU 

Ravelet et a/.ll2010HPaoletti fc Lathropl201ll) and numerical data ()Pirro fc Quadriol200 
Bils on fc Bremhorstll2007l:lDonej 120071 12008) are shown. The dashed lines are Esser and 
Grossmann's (1996) estimate for the onset of turbulence with T) — 0.71. The dark shaded 
area indicates the data points in the well studie d small Reynold s number regime of pattern 
forma ti on and spatial temporal cha os (see e.g. Andereck et al. ( 19861 ); Pfister fc Rehberg 
( 198lh : ICross fc Hohenberd l| 19931) 1. The light gray area is the region shown in Figure 
2"al covered by the present DNS. 



conditions also in azimuthal direction and an aspect ratio r = 2 in axial direction rather 
than r = 2ir as we do. Also they find a maximum in the angular velocity transport 
for moderate cou n ter-rotation a = —u /uji = 0.4, similar as found in the experiments 
bv Ivan Gils et at ( 20116 . 2012 ) and in the present numerical simulations. So the result 
seems to be very robust and does at least not strongly depend on Ta, T and other details 
of the flow. iBrauckmann fc Eckhardtl ([20121 ) also offer an analysis of the PDFs of the 
local angular v elocity fluxes in the d ifferent regime, similarly as has been done in the 
experiments bv lvan Gils et al. ( 2012 ). 



2. Numerical method 

The Taylor- Couette flow was simulated by solving the Navier-Stokes equations in a 
rotating frame, which was chosen to rotate with = uj e z . This way the boundary 
conditions are simplified: at the inner cylinder the new boundary condition is ug(r = 
r i) — Tifai — cj ), while at the outer cylinder we have a stationary wall ug(r = r Q ) = 0. 
We can choose the characteristic velocity U = r^uji ~ oj \ and the characteristic length 
scale d to non-dimensionalize the equations and boundary conditions. The characteristic 
velocity U can be written as 

U = (u/d)-[8r, 2 /(l + r ] f]-Ta 1 ^. (2.1) 

Up to a geometric factor, which is 0.810 for our choice of rj, the characteristic velocity 
U is the Ta 1 / 2 -fold of the rotation, independent of the molecular velocity vjd. The non- 
dimensional variables will be labeled with a hat. In this notation, the non-dimensional 
inner cylinder velocity boundary condition simplifies to: ug(r = ri) = (cj s ; — u )/\uji — w \. 
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As LJi — uj is positive throughout the range covered in this work, in our coordinate system 
the flow geometry is simplified to a pure inner cylinder rotation with the boundary 
condition Ug{ri) = 1. The effect of the outer cylinder is felt as a Coriolis force in this 
rotating frame. 

The resulting Navier-Stokes equations in the rotating frame are now 

—t- it • V"U = -Vp + (^] ' V 2 m + Ro~ 1 e z x u , (2.2) 



dt \Ta 
where the Rossby number is defined as 



Ro= \^Ar \l + a\ r, 

2uj d a 2(1 - T}) v ' 



and f(r)) as 



877 

Equation (|2.2I) is in analogy to the Navier-Stokes equation for a rotating Rayleigh- 



m - . (2.4) 



Benard system, 



„ +«■ Vu = -V/3 + — V 2 tt + 6»e z - i?o _1 e z x m , (2.5) 

with the main difference that the Rossby number's sign (carried by uj in eq. (|2.3[1 ) is 
relevant in Taylor-Couette flow. As long as the transport of angular momentum takes 
place from the inner to the outer cylinder, i.e. u>i > w D , Ro is always negative for counter- 
rotating cylinders and always positive for co-rotating cylinders. Therefore the sign of Ro 
affects the flow physics, as it indicates the direction of rotation of the outer cylinder. 

These equations were solved using a finite difference solver in cylindrical coordinates. 
The domain was taken to be periodic in the axial direction. Coordinates were distributed 
uniformly in the axial and azimuthal direction. In the radial direction, Chebyschev clus- 
tering was used to cluster points near both walls. For spatial discretization, a second 
order scheme was used. Time integration was performed fractionally, using a third order 
implicit Runge-Kut t a met hod. More details of the numerical method can be found in 
Verzicco fc Orlandl (1996). Large scale parallelization was done using a combination of 



MPI and OpenMP directives. 

In order to quantify the flow, it is useful to continue with the normalized radius f = 
(r — f*i)/(r — r,-) and the normalized height z — zj \r — r,j). As an aid to quantification, 
we define the time- and azimuthally-averaged velocity field as: 



&(r,z) = (u{9,r,z,t)) e ,t . (2.6) 

This time- and ^-independent velocity is used to quantify the large time scale circulation 
through the wind Reynolds number: 

Re w = with U w = U(u r + u z )\%. (2.7) 

As mentioned in eq. (|2.1I) . U cx Rei — rjRe scales as Ta 1 / 2 ; the non-dimensinal transverse 
velocity fluctuations may or may not lead to corrections of this basic scaling. 

The convective dissipation per unit mass can be calculated either from its definition 
as a volume average of the local dissipation rate for an incompressible fluid 
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Ret Nu w 


(present study) Nuu, 


fFasel & Booz 1984) Nu„, 


(Pirro & Quadrio 2008) 


60 


1.0005 


1.0000 


1.0000 


68 


1.0006 


1.0000 


1.0000 


70 


1.0235 


1.0237 


1.0238 


75 


1.0835 


1.0833 


1.0834 


80 


1.1375 


1.1371 


1.1372 


Table 2: N Ulu 


for low Rei number and Re — 0, T = 2, r\ = 


^ and rotational symmetry 


of order 4. 









e u = V ({Vuf) vt (2.8) 
or from the global balance (EGL 2007): 

e u - e ufi = -^a~ 2 Ta{Nu u - 1) , (2.9) 

where e Ut o is the volume averaged dissipation rate in the basic, azimuthally symmetric 
laminar flow, cf. eq. (I1.16[) . 

In order to validate the code we have calculated e u from both (|2.8p and (|2.9[) and 
checked for sufficient agreement. We define the quantity A e measuring the relative dif- 
ference 



/ i d- i a- 2 Ta{Nu LJ - 1) + e ufi - i/((Vu) 2 ) 



^ . (2.10) 



^((V«) 2 ) y , t 

A e is equal to analytically, but will deviate when calculated numerically. 

The strictest requirement for numerical convergence was that the radial dependence of 
Nuufr) had to be less than 1%. This is a much harder condition to satisfy than torque 
equality at both inner and outer cylinders, which is satisfied if the Nu^ at both cylinders 
is equal. Indeed, in many cases the torques were equal to within 0.01% but Nu^r) was 
not constant within 1%, which meant either a higher resolution had to be chosen or that 
the simulation had to be run for longer time. The time-average of the energy dissipation 
calculated locally (equation 12. 8p was also checked to converge within 1%; see section [ 
for more details. 



3. Code validation 

3.1. Validation against other codes at low Reynolds number 

Firs t of all, the code was va lidated against numerical results from iFasel fc Boo3 (Il984l) 
and Pirro fe Quadrio (2008) . This comparison was done through Nu^ measurements at 
low Rei numbers, in the range between 60 and 80. Only a quarter of the TC system was 
used, assuming a rotational symmetry of order four. The aspect ratio T was taken as 
two, the radius ratio r] as 0.5. These geometrical parameters are different than the ones 
used in the rest of the paper, but they are used here for validation. The resolution of 
the sim ulation (Ne x N r x N%) was taken as 32x64x64, the same as in Pirro fc Quadrio! 
(2008). The results can be seen in Table[2j The values show a match up to three significant 
figures, or sometimes even higher. 

For the two smallest Reynolds numbers, both references obtain the same result, while 
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Fi gure 4: Nu^, versus Ta for r\ — 5/7. Experiments (circles) and numerics agree in shape, 
but there is a slight shift between the data, which we attribute to the different boundary 
conditions in lateral direction. 

we obtain a slightly different result. This probably comes from the fact that they measure 
the torque directly at the inner cylinder, which we then convert to Nu u for comparison, 
while we measure Nu u by taking an average value of J{r) and converting this to a value 
for the torque and thus Nu^. The difference between these approaches is probably the 
origin of the discrepancy. However, as it is very small (below 0.1%) it is not worrying. 



The code was als o validated by comparing responses obtained at higher Taylor numbers 
versus data from lLewis fc Swinnevl (J1999) . This was done through the Nusselt number 
for pure inner cylinder rotation (a = 0) at fixed rj = 5/7 and T = lix. The overlap 
between the simulations and experimental data can clearly be seen in the higher Taylor 
number range, which we have achieved with the numerics. The shift of about 5% might 
be attributed to the difference in both aspect ratio and boundary conditions at the top 
and bottom because of the vertical confinement in the experiment. As we also have an 
overlap at the lower Taylor range with other numerical simulations as shown in Section 
13.11 we feel sufficiently confident to proceed with our code. 

3.3. Resolution tests 

The grid spatial resolution was tested at a = by calculating Nu^ beyond onset of Taylor 
vortices (Nu^ > 1) and A e from eq. (|2.10l) . which analytically is equal to 0. Both was done 
for different grid resolutions with increasing Taylor number. For all these simulations we 
took r = 2tt, r\ = 5/7 and a = 0. The results are shown in Table [3] 

Spatial convergence required more grid points than initially expected as satisfying 
the torque balance is a necessary but not a sufficient condition for grid independence. 
The top two graphs in figure [5] show a plot of the radial dependence of Nu^ at Ta = 
2.5 • 10 5 for an underresolved case (100x50x50, Nu u = 2.70845), reasonably resolved case 
(200x100x100, Nu u = 2.76208) and extremely well resolved reference case (300x150x150, 
Nu u = 2.77855). Nu u should not be a function of the radius as mentioned previously, 
but it does show some dependence. For the underresolved case we can see that the torque 
balance is satisfied very well (0.06%), even if other criteria are not satisfied as the peak- 
to-peak variation of Nu u is approximately 1% and the relative error in comparison to 
the reference case is 2.5%. The graph also shows that taking the value of Nuuj at one of 
the cylinders gives a higher result than taking the radial mean. 



3.2. Comparison with experiment 
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Ta 


N e x N r x N z 


Nu u 


-A e 


Case 


4 ■ 


10 4 


128x64x64 


1.86927 


0.0159 


Resolved 


4- 


10 4 


256x128x128 


1.85562 


0.0074 


Error Reference 


1 • 


10 5 


160x80x80 


2.40536 


0.0215 


Resolved 


1 • 


10 5 


256x128x128 


2.40216 


0.0322 


Error Reference 


2.5 


■ 10 5 


100x50x50 


2.70845 


0.0392 


Underresolved 


2.5 


■ 10 5 


200x100x100 


2.76208 


0.0102 


Resolved 


2.5 


• 10 5 


300x150x150 


2.77855 


0.0062 


Error Reference 


7.5 


• 10 5 


256x128x128 


3.49816 


0.0147 


Resolved 


7.5 


• 10 5 


384x192x192 


3.51268 


0.0056 


Error Reference 


2 • 


10 6 


192x96x96 


4.83540 


0.0949 


Underresolved 


2 ■ 


10 6 


256x128x128 


4.46000 


0.0174 


Resolved 


2 ■ 


10 8 


384x192x192 


4.47765 


0.0065 


Error Reference 


4- 


10 8 


300x144x144 


5.42553 


0.0216 


Resolved 


4 ■ 


10 8 


432x216x216 


5.37264 


0.0063 


Error Reference 


1 • 


10 7 


384x192x192 


6.42160 


0.0168 


Resolved 


1 • 


10 7 


641x321x321 


6.34068 


0.0078 


Error Reference 


2.5 


•10 7 


641x321x321 


7.46617 


0.0161 


Resolved 


5- 


10 7 


800x400x400 


8.76601 


0.0166 


Resolved 


1 ■ 


10 s 


1024x500x512 


10.4360 


0.0170 


Resolved 



Table 3: Resolution tests for r\ = 0.714 and T = 2ir. The first column displays the 
Taylor number, the second column displays the resolution employed, the third column the 
calculated Nu u , the fourth column the relative discrepancy A e between the two different 
ways of calculating the energy dissipation, and the last column the 'case'. The resolution 
is sufficient for all cases, as the variations are small. A e turns out to be negative; thus 
the code gives for the dissipation rate a larger value for the derivatives-squared-based 
definition than for the iVit^-based balance expression. 



The bottom two panels in figure [5] show the same plots for Ta = 2 ■ 10 6 and the three 
cases: underresolved (192x96x96, Nu^ — 4.8354), reasonably resolved (256x128x128, 
Nuu = 4.4600) and extremely well resolved reference (384x192x192, Nu^ = 4.4776). 
For this Taylor number, the underresolved case shows a smaller deviation of Nuuj from 
the mean value and the torque difference when compared to the lower Taylor number. 
However, the discrepancy in the mean value of Nu^ between the underresolved and the 
reference case is much larger (7.9%). For this Ta the value of Nu^ at the cylinder walls 
is larger than the average value of Nu^, too. 

If we look at the Nu u (r) profiles at given Taylor numbers Ta, they show similar radial 
dependences, whose magnitudes decrease with increasing resolution. However, the shape 
of this dependence is different for both Taylor numbers. The peaks of 7Vw w (r) are located 
close to the boundaries, indicating that they are probably produced by some boundary 
layer features and are not a systematic bias of our solver. 

Dissipation should be equal, irrespective of the way in which it is calculated, directl 
from its definition or indirectly via the Nusselt number balance. IStevens et all (|201 
also mentioned the importance of the corresponding equality in Rayleigh-Benard flow 
for low values of Pr as a way to ensure that the flow field is sufficiently resolved and 
that the gradients are captured adequately. Underresolving a flow in Taylor- Couette will 
result in a value of A e which is too large in magnitude. This can be seen in Table |3] 
for the underresolved simulations at Ta = 2.5 • 10 5 and Ta = 2 • 10 6 . We consider the 
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Figure 5: Radial dependence of Nu^f) for three different grid resolutions (see legends). 
The top two figures are for Ta — 2.5 • 10 5 and the bottom two are for Ta — 2 ■ 10 6 . 
The figures on the left show the absolute values, an error bar indicating a 1% error for 
reference. The resolved cases lie within this error bar. The figures on the right show the 
curves normalized by their average value to compare the radial fluctuations of Nu^. 



convergence of A e towards zero as a neccesary but not a necessarily sufficient way to 
ensure grid convergence. 



3.4. Dependence on initial conditions 

For the lower Taylor numbers Ta the flow was started from rest (it = 0). The Taylor 
vortices start forming within a couple of revolutions. After enough time, a steady state 
with three pairs of Taylor vortices was reached. However, the simulations can also be 
started from non-resting conditions. Depending on these conditions a different number 
of vortex pairs can arise. This has a strong influence on both the global and the local 
response of the system. Once the vortices have formed, they are persistent in time during 
the simulation. Therefore, it is possible to bias a simulation through the initial conditions 
to have a higher or lower amount of vortex pairs, which results in a different response. 

Although the importance of these coherent structures gets smaller and smaller with 
increasing Ta (Section [5]), at lower Ta a certain number of vortex pairs must be fixed to 
determine the response. The three vortex pair state has been chosen to be the base state, 
to keep the aspect ratio of the vortices as close to 1 as possible (2N pa i rs /T = 0.96). The 
effect of having 3 or 4 vortex pairs on the response for selected To is shown in Figure 
[6] It is important to note that this effect is different from the effect caused by neutral 
surface stabilization, which occurs when the vortices cannot penetrate the whole flow, 
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Figure 6: Dependence of Nu u — 1 vs Ta on the number of vortices. Although the values 
of Nuoj are different, the scaling behavior with Ta is the same. The error bar indicates 
a 5% difference. 



and the number of vortices is changed as a result. That will be featured in more detail 
in Section IB~5l 



4. Global response 

In this section, the global response of the Taylor- Couette system is shown across the 
parameter space. First the onset of Taylor vortices is analysed. Then the scaling laws are 
revealed for pure inner cylinder rotation. Finally, the effect of the outer cylinder rotation 
on the scaling laws is investigated and an optimum of Nu^ as a f unction of a for given T a 
is found, as has been reported for large Ta from experiment, cf. Ivan Gils et al. ( 2012 ). 



4.1. Transport Nu u and wind Re w for pure inner cylinder rotation 

The global response of the system is quantified through Nu u and Re w . These two quanti- 
ties measure two different flow responses. Nu u quantifies the transport of angular velocity 
and Re w the "wind" , i. e. the additional velocity on top of the azimuthal flow. For the 
purely laminar-azimuthal flow Nu u = 1 by definition, and Re w = as this laminar flow 
only has an azimuthal velocity component. 

First of all we analyse how the onset of Taylor vortices is reflected in the global response 
quantities Nu u and Re w . Nu u — 1 is is the additional transport of angular velocity on top 
of the laminar transport and the wind is the fluid motion on top of the purely laminar- 
azimuthal flow. Figure [7] shows the numerically calculated Nu u — 1 and Re w as functions 
of Ta close to onset of the Taylor-vortex state. The critical Taylor number (Ta c ) for 
the onset of Taylor vortices is calculated to be around 1020 for our value of rj. This 
DN S value can be compared with Ta c as obtained from the analytical approximation 
of Esser fc Grossmannl (1996), which is 1038 for the present rj. The agreement is within 



1.6%. Later on we shall use these analytically calculated onset Taylor numbers. 

After the Taylor vortices have appeared in the system, they are the dominating feature 
of the flow for several decades of Ta. The top two panels in figure [8] show the response of 
the system with increasing Taylor number in the case of resting outer cylinder and pure 
inner cylinder rotation. We plot Nu^ vs. Ta-Ta c rather than vs. Ta as it then shows a 
better scaling for the points at low Ta. 
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Figure 7: System responses Nu^ — 1 (left) and Re w (right) as a function of Ta for pure 
inner cylinder rotation near the onset of Taylor vortices. Onset in the present DNS occurs 
at Ta fa 1020. 



There seems to be a clear change in the scaling law of Nu^ versus Ta, but not so 
in the scaling law for the wind Re w as a function of Ta. This change occurs between 
Ta = 2 • 10 6 and To = 4 • 10 6 and has been seen in other numerical simulations too 
(|Coughlin fc Marcuslll996h . Nu u - 1 scales as {Ta - Ta c ) 034 for Ta < 2 x 10 6 and as 
Nu u - 1 ~ (Ta - Ta c ) - 21 for Ta > 2 x 10 6 . We attribute the change in scaling to the 
changes in the coherent flow structures that affect the heat transport but not the global 
wind amplitude. As will be discussed later in detail, we expect coherent flow structures 
to lose importance for increasing Ta, see Section [5] We note already here that although 
the loss of influence of coherent structures in RB convection (for Pr = 1) and in TC 
flow sets in at similar values of Ra and Ta respectively, i.e. around 10 7 ( Sugivama et all 
2007), there is a large difference in the shear Reynolds numbers Re s of the boundary 
layers in these two systems. This will be discussed in the following Section. 

Re w measures the amplitude (strength) of the Taylor vortices, which persist at long 
time scales. The nondimensional characteristic speed for these vortices remains approxi- 
mately constant with Ta, viz. about 5-6% of the inner cylinder velocity Uj, throughout 
the whole Taylor number range considered, and that is why we see a direct scaling law 
of Re w ~U ~ (Ta- Ta c )?, cf. ea.(l2~11. 

The mutual functional dependence of the two responses Nu u and Re w is presented in 
the bottom left panel of figure [5] As expected from Figures I8aj|8bl the relation between 
Re w and Nu^ also shows the change in the scaling. We interprete this as follows. Be- 
fore the change, mainly the Taylor vortices are responsible for the additional transport. 
Beyond the change, some short time scale fluctuations appear, indicating other coherent 
structures, which disrupt the flow and finally become its dominating features, while the 
Taylor vortices lose importance. In order to see these time scales, we show the temporal 
dependence of N ii u (t) in Figure [8fJ 

The Nusselt number shows almost no time dependence for lower Taylor numbers. But 
it shows two different time scales at Ta — 4 • 10 6 . The short time scale gains much more 
importance for the highest Taylor numbers, causing fluctuations of about 10%. Figure [9] 
shows Nu(t) for two Ta numbers between which the transition to time dependence takes 
place. For Ta = 2.75- 10 6 the global response is still stable to perturbations and the time 
dependence dissapears with time. However, for Ta = 3 • 10 6 , these perturbations do not 
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Figure 8: The top two figures show the system response as a function of Ta with best 
fit lines for pure inner cylinder rotation, Nuu, — 1 on the left and Re w on the right. The 
middle plots show the date of the two top plots compensated by the scaling factor to 
test the quality of the scaling behavior. The bottom left figure presents the functional 
relation between the two responses Nu^ — 1 and Re Wl and the bottom right figure shows 
the temporal dependence Nu^it) for three different Ta. The time dependence can be seen 
to set in between Ta = 2-10 6 and Ta = 4-10 6 , just where we s ee the change in the effective 
scaling in the two left figures. The analytical approximation of Esser fc Grossmann ( 19961 ) 
is used, i.e. Ta c — 1038 for the present rj. 



decay, and a small constant oscillation remains. This oscillation becomes larger and more 
chaotic with increasing Ta as seen previously in Figure [8f| 
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Figure 9: Nu u (t) behavior with time at Ta = 2.75 ■ 10 6 (left) and Ta = 3 • 10 6 (right), 
Nuu,(t) can be seen to be stable to perturbations at the lower Ta. At the higher Ta, a 
small oscillatory behaviour remains. At Ta = 4 • 10 6 and beyond the time dependence is 
strongly chaotic, cf. Figure l8fl 



4.2. The effect of outer cylinder rotation and optimal transport 

In this section the effect of the outer cylinder rotation on the global responses Nu^ and 
Re w will be studied. This effect is felt by the flow as a Coriolis force (Equation \2.2\i , so 
plots in this section will be done versus i?o _1 oc — a/|l + a\ with a = —uj /uji. 

Ro -i = ™*L = -2^^ = 2 l -\ u ° , . (4.1) 
U r) \l + a\ r\ \u)i—(jj \ 

The inverse Rossby number i?o _1 runs from 277/(1 + 77) to —1 if u) runs from rfuji, the 
inviscid Rayleigh-line in the first quadrant of the (Re , i?ei)-plane, to —00. It is useful to 
note that for a given Ta constant i?o _1 means constant a and vice versa. Thus a seems 
the preferrable notion as it is more direct; i?o _1 will be used only, when it provides a clear 
advantage in visualization or later in the paper when we will trace badk the occurrence 
of the maximum to the Navier- Stokes equation. 

Figure [TU] shows the complete results for Nu u as a function of i?o _1 and Ta. In order 
to better quantify the results from Figure [TUl cross sectional cuts are taken. By taking 
cross sections of constant a, scaling laws can also be discovered for non-zero values of 
a, i.e. for co- and counter-rotation ui a 7^ 0. This is shown in Figure [TT] for five different 
values of a, two under co-rotation lo > 0, two under counter-rotation uj < 0, and as 
reference case uj d = 0. 

For counter-rotating cylinders and Taylor numbers prior to the change in scaling hap- 
pens, a universal scaling of approximately Nu u — 1 ~ (Ta — Ta c )°- 34 is seen. However, 
the change in scaling and its exponent happens earlier in Ta for a = 0, while the scaling 
prevails to larger Ta for the other a (0.2 and 0.4). The scaling is different for co-rotating 
cylinders, and is approximately Nu^ — 1 ~ (Ta — Ta c )°- 27 . Ta c = 1038 has been sub- 
tracted as done previously so that the scaling is not lost for the first points. 

The time independence of JVm u is broken for much smaller Ta, if the outer cylinder is 
rotating. For both co- and counter-rotating cylinders time dependence can be noted to 
set in at Ta as low as 10 5 . Also, the scaling of Nu^ with Ta is maintained throughout a 
much larger range of the Taylor number. Therefore, the breakdown of time independence 
cannot be associated anymore to the change in scaling, as one could conclude when only 
considering pure inner cylinder rotation, where the loss of time-independence and change 
in scaling happened at about the same Ta. 
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Cross sections of constant Ta are shown in Figure Q21 Nu u — 1 = indicates points 
for which the flow is purely laminar-azimuthal. For co-rotating cylinders, the maximum 
value of -Ro -1 reaches the inviscid Rayleigh stability line Ro]^ = 0.833 for even the lowest 
values of Ta. On the other hand, the minimum Ro^ 1 which destabilizes the laminar state 
can be seen to decrease (become more negative) with increasing Ta. 

In order to better visualize the results it is useful to define a normalized Nusselt number 
as Nuu = {Nu u {Ro^ 1 ) — l)/(Nu UJ (Ro^ 1 = 0) — 1), which allows easier visualization of 
the dependence of Nu u on i?o _1 across the Ta range of interest. The numerator of Nu u 
goes to zero, if Ro~ l becomes too large or too small, i. e. reaches the stability lines 
(where Nu u becomes 1, since the flow is laminar-azimuthal in the stable ranges), while 
the denominator is always larger than zero, as long as Ta > Ta c . 

If Ta is large enough the shape of the graph resembles two straight lines from the 
maximum value of Nu^ to Ro~^ a and i?o~ in . These straight lines have already been 



seen w hen plotting Nuu, versus a slightly different version of i?o _1 in lPaoletti &: Lathrop 
( 201 lh and are the reason we chose to plot Nuu versus Ro^ 1 in this section. 

For smaller Ta an optimum transport for co-rotation can be seen, i. e., for positive 
values of Ro^ 1 . This holds for Ta less than the discussed change of the scaling behavior of 
Nuu versus Ta. For the two lowest values of Ta the deviation of i?o _1 beyond may still 
be within numerical uncertainties. However, Ro~^ t is definitely positive for Ta = 7.5 ■ 10 5 
and Ta — 2 ■ 10 6 and seems to fit with the piecewise linear shape of the graphs. At 
around Ta = 4 • 10 6 , which is beyond the change in scaling, the maximum Nu^ begins to 
drift to negative Ro -1 , i. e., towards counter-rotation. This can be seen in Figure IT3"1 In 
fig. [14] we plotted the positions of the optimum transport in the (i?e Q , Re{) phase space. 
Clearly, the curve does not have equal distance to the two instability branc hes of the 



Esser-Grossmann approximation, as was speculated in Ivan Gils et all (|2012h . Another 



feature of the drift of Ro Q } t is the following: While the curve of Nu u has a prominent 
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Figure 11: Nu^ — 1 versus Ta — Ta c for three values of co-rotating a (left) and pure inner 
cylinder and two values of counter-rotating a (right). Compensated plots are shown below. 
Numerical errors are less than 1%. Ta c depends on a and is determined respectively from 
the analytical approximation ( Esser fc Grossmannlll996l) ; for a=0 it is Ta c = 1038. 




Figure 12: Nu^ — 1 as a function of i?o _1 across the available Ta range. The advantage of 
plotting versus i?o _1 can be seen as some of the graphs show piecewise linear behaviour 
and the value of Ro^ 1 at the inviscid Rayleigh stability limit a = —rj 2 (denoted as Ro^) 
is Ro^ = 0.833, appears for all Ta. 





Figure 14: Location of the optimal transport (black dots connected by thick black line) 
in the (Re ,Rei) phas e space of figure Ek- The s olid blue line represents the onset of 
instability according to Ess er fc Grossman nl (|l996h and the thin dashed line is the line of 
equal distance to the left and right branch of the Esser-Grossmann instability line. The 
solid line is the bisector of the Rayleigh instability line (a — —rj 2 ) and the line of pure 
outer cylinder rotation (a — oo). 



peak at Ro^ for values of Ta of around 10 6 — 10 7 , this turns into a plateau for the 
highest value of Ta and Ro~^ t becomes hard to identify. For the highest value of Ta, the 
lower border Ro ~)„ already is beyond our parameter range of negative Ro~ x . 

Experiments ( van Gils et al. 20116f ) have found an optimum transport a op t ~ 0.33, 
corresponding to Ro~^ t ~ —0.20 for Taylor numbers of the order of 10 12 . Thus the 
position of the maximum shifts for higher Taylor numbers. 

Figures [12] and [13] show some anomalous jumps in the graph around Ro^ 1 w —0.2 
which corresponds to a sa 0.3. These are caused by different vortical states as mentioned 
in section 13.41 These may be present even if the simulations start from the same initial 
conditions for different values of Ro~ l and Ta. If the number of vortices is higher, the 
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vortices become stronger, and the value of Re w which measures their strength also be- 
comes higher. Since Nu u is monotonously related with Re w , it also increases. We will 
further analyse this multi- vortex state in Section 16.51 



5. Characterization of the flow state 



In this section we will analyse two characteristic Taylor number ranges in TC flow. 
The first, lower one, is the range in which the importance of coherent flow structures is 
lost, since these have become too small in size. In section 2] we have observed a change in 
the scaling law for the heat flux from Nu u ~ (Ta - Ta c ) 34 to N Ucj ~ (Ta - Ta c ) a21 . 
Although the Taylor number for this change coincides with the onset of time dependence 
for pure inner cylinder rotation, when adding outer cylinder rotation the onset of the 
time dependence is much earlier, and a transition in the scaling laws cannot even be 
seen. 

Therefore, the onset of time dependence cannot be linked satisfactorily to the change 
in the -/Vu^-scaling. Another way of explaining this change is by estimating when the 
spatially coherent flow structures loose influence because their size becomes too small. 
We do this b y defining an average , global coherence length in terms of the Kolmogorov 
length scale ([Sugivama et alluOOlv resulting from the volume averaged dissipation rate: 



4 = 10 VK = 10(z/ 3 /e«) 1/4 = l0d{a- 2 Ta{N UbJ - 1) + e u , )- 1/4 . (5.1) 

where eq. (|2.9p has been used for the second equality. We compare the global coherence 
length £ c with the gap width d or, equivalently with the extension of the remnants of 
the Taylor vortices, whose length can also be estimated as d, since they tend to have a 
square aspect ratio. 

Figure [15] shows the variation of l c jd with increasing Taylor number. The loss of 
importance of coherent structures happens in the range where £ c /d is between 0.1 — 0.5, 
corresponding to Ta rs 10 6 — 10 7 . It is just within this Ta range, where the change in 
the Nu u scaling law occurs. The graph is consistent with that change taking place at 
approximately the same Ta for different values of a, which is what is seen in Figure ITTl 

A second characteristic Taylor number is connected with the shear instability of the 
boundary layers (BL). Here the laminar Prandtl type BLs become turbulent. Beyond that 
T a the flow is fully turbulent throughout and this state is known as the ultimate state, 
cf.[G rossmann fc Lohsel (|201lh . This happens if the BL shear Reynolds n umber becomes 



Re s > 250 — 420 and Re s is calculated from the shear velocity U s as ([van Gils et 



2012) 



Re s 



U S S 



apB 



VRe 



Re v 



(5.2) 



The empirical constant aps is taken as 2.3 for Prandtl-Blasius type boundary layer in 
TC flow. 

Fo r RB convect i on th i s transition is expected at Ra w 10 14 (jGrossmann fc Lohse 
20011: lAhlers et all 120091 ; IGrossmann fc Lohsel 12011 ). while figure [T5bl shows that the 
transition in TC is expected for 10 s < Ta < 10 9 . Recently, experim ents have confirmed 
the ultimate scaling both for Nu, 
TC flow N Uu] ~ Ta 38 and Re w - 



and Re w . iHuisman et al. ( 20121 ) have shown that in 
Ta 50 when Ta > 10 9 . A confirmatio n of the analogy 
between RB and TC is obtained by the high Ra number experiments bv lHe et al\ ( 201 lh 
as they measured that Nu ~ Ra 3S and Re w ~ Ra 50 fo r Ra > 5 x 10 14 . These measu red 
scaling exponents agree exactly with the predictions by IGrossmann fc Lohsel ( 2011 ). In 
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Figure 15: The graph on the left shows £ c /d versus Ta for i] — 5/7 and three values of a 
as shown in the inset. The blue region indicates the range of decreasing importance of the 
coherence structures. - The graph on the right shows the shear Reynolds number Re s of 
the boundary layers versus Ta. In the blue shaded region we expect a shear instability 
of the boundary layers, in which the Prandtl laminar bound ary layers become turbulent . 
The TC system then is in the so-called 'ultimate' regime, cf. Grossmann fc Lohsel ( 20 1 lh . 



contrast to the experiments of Ivan Gils et al\ (|20116l . 120121 ). in our present numerical 
simulations the ultimate state is not yet achieved, as clearly seen from fig. I15bl 



6. Local results 

This section contains the analysis of local results. For convenience we skip in this 
section the "hat" on the dimensionless flow field variables, but still understanding them 
as being dimensionless. The angular velocity profiles are shown and the ratio of the 
BL thickness is calculated and compared with the theory of EGL (2007). The angular 
velocity profiles reflect the interplay of bulk and boundary layers and that of the mean 
flow and added perturbations. The importance of convective versus diffusive transport is 
quantified through the bulk slope of the angular velocity profile, and again we will find 
a maximum as function of a, which we will connect with the maximum in the angular 
velocity transport Nu u . 

6.1. Local coherence length 

Figure[l6]shows the local coherence length calculated from the local dissipation in analogy 
to eq. (I5.1[) . This figure adds details on where we expect the Taylor vortices to break 
down. At low Taylor number, the local coherence length is smaller than 0.1 only very 
near to the wall, where the highest local dissipation takes place. With increasing Taylor 
number, the highest local dissipation still is near the wall, but the dissipation rate is 
large enough in the whole domain to break up the dominance of the coherent structures, 
even if they do not fully disappear but become small enough. 

6.2. Angular velocity profiles 

The angular velocity w is the transported quantity in Taylor- Couette flow. Analysing 
the dependence of the io(r) profiles on the driving parameters Ta and a seems useful 
to understand how transport takes place in the flow. uj{r) profiles are shown in Figures 
1171 and 1181 Beyond the breakdown of the laminar, purely azimuthal flow, three distinct 
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Figure 16: The two plots show the areas in red where £ c /d < 0.1 for pure inner cylinder 
rotation and Taylor numbers of (a) Ta = 2 • 10 6 and (b) Ta = 1 • 10 7 . In these areas the 
breakup of the coherent structures is likely to occur. The arrows indicate u in magnitude 
and direction. 



regions in the gap appear. These are the inner and outer boundary layers (BL), in which 
the transport mechanism is dominantly diffusive, and a flatter bulk zone, in which the 
transport mechanism is dominantly convective, see Figure I20bl for a sketch, in which we 
approximate the profile of the mean azimuthal velocity (ug) z by three straight lines, one 
for each boundary layer and one for the bulk. For the boundary layers, we calculate the 
slope of the lines by fitting (least mean square) a line through the first three computa- 
tional grid points. For the bulk, we first force the line to pass through the inflection point 
of the profile (the nearest grid point). Then, its slope is taken from a least mean square 
fit using two grid points on either side of this inflection point. The respective boundary 
layer line will cross with the bulk line at a point which then defines the thickness of that 
boundary layer. 

In the next two subsections we will discuss the BL and bulk regimes separately. 

6.3. Angular velocity profiles and resulting boundary layer thicknesses 

With increasing Ta, in order to accommodate the increasing angular velocity transport, 
the boundary layers become thinner and the w-slopes steeper. What is striking is the 
strong asymmetry between the inner BL and the outer BL, which is much thicker. Figure 
[T9l shows the ratio between the outer and the inner boundary layer thicknesses versus 
rotation ratio a for four Taylor numbers. This asymmetry is a consequence of the exact 
relation d r (ui}\ = r] 3 d r (uj}\i, obtained from the r-independence of J", cf. eq. (|1.4I ): The 
slope at the inner cylinder is a factor of ?/~ 3 larger than at the outer one and thus the 
outer boundary layer is much more extended than the inner one. Since for the present 
Ta-range the shear Reynolds number Re s is still below the threshold value range for 
the transition to turbulence in the boundary layers (see section [5]) , we can compare the 
numerically obtained boundary thickness ratio with that one obtained by EGL (2007), 
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Figure 17: The t,8 and also z-averaged angular velocity (uj) z versus f for four Taylor 
numbers. The boundary layers of the U) profiles become thinner around slight counter- 
rotation than they are for higher values of a as well as for co-rotation as an indication 
of increased transport. For strong co-rotation a = —0.4 as well as high a, i. e., strong 
counter-rotation, at low Ta there is not yet a flat bulk zone since the flow is not yet 
turbulent enough. 



which had been derived in the spirit of the Prandtl-Blasius (i.e., laminar- type) boundary 
layer theory, namely 

K = - 3 \<»o-0>Mk\ (61) 

\Ui-Ubuik\ 

Here the value of LUbuik is a characteristic bulk angular velocity chosen to be the value 
of lu at the inflection point of the w-profile (see figure [20]) . It is calculated from the 
numerical simulations. The result for the BL thickness ratio A°/A^ is shown in figure 
[Hfl The agreement with the numerically obtained ratio is satisfactory for the counter- 
rotating a-cases, getting even better with increasing Ta. This is because the estimate is 
based on a flat profile in the bulk, and indeed the profile becomes flatter with increasing 
Ta. For co-rotation, formula (|6.ip apparently fails. This had to be expected, because the 
approximation of the profile of (ug) z by three straight lines, which was assumed in the 
derivation of (|6.1j) . is then no longer appropriate. 



6.4. Angular velocity profiles in the bulk 

We now come back to the mean profiles in the bulk. As can be clearly seen from comparing 
figs. [17] and [HI both the mean angular velocity and its slope are controlled by a (or i?o _1 ) 
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Figure 18: (cD) 2 versus f for four values of a and increasing Ta. The boundary layers 
for the uj transport become thinner with increasing Ta indicating increased angular 
velocity transport. The boundary layers get steeper with increasing Ta and the bulk 
region becomes more extended. For low values of Ta the bulk region is rather small. The 
w(f) profile in the center remains approximately unchanged with Ta. 

rather than by Ta. This behavior can be understood from equation (I2.2[) : The outer 
cylinder rotation is reflected in that equation as a Coriolis force. This force is present in 
the whole domain, while Ta controls the strength of the viscous term, which is dominant 
in the boundary layer. Therefore the profile is controlled by the Coriolis force, i.e., i?o _1 
or a, and not by Ta. 

To further quantify this, the gradient of (ui) z is calculated. This is done by fitting a 
straight line to (u) z (f) at the point of the profile's inflection, numerically using the grid 
points around it. An example of how this is done can be seen in Figure [2Ual The results 
for the profile slopes in the bulk as functions of Ta and a or i?o _1 are shown in Figure 

HO 

The graphs collapse on each other for co- rotation (a < 0), which is what we expected 
from Figure IT51 and our previous analysis. An almost linear relationship between i?o _1 
and the bulk slope of w(r) is found. If pure inner cylinder rotation (a = 0) is approached, 
the graphs for the various Ta start to differentiate and reach a plateau. The absolute 
value of the slope decreases with increasing Ta. This is due to the increasing impor- 
tance of convection at higher Ta. Note however that also this counter-rotating case, for 
large enough Ta the center slopes again lose their Ta dependence, i.e., are again mainly 
controlled by i?o _1 and thus the Coriolis term. 
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Figure 19: A°/A^ versus a for four values of Ta. The agreement between theory and 
simulation is better for counter-rotation and with increasing Ta, but does not match for 
co-rotation. A dashed line indicates the optimal rotation ratio a op t- 




Figure 20: Example of the two fitting procedures for Ta=2 ■ 10 6 and a=— 0.2. The left 
panel shows a linear fit to the bulk part of the angular velocity (w) 2 . The right panel 
shows in addition the fit to (ug) z for its boundary parts, which will be used in Section 
16.31 Both bulk fits are done at the inflection point, but for different variables (ui or ug), 
which gives a slight difference. 
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Figure 21: The average slope of the angular velocity d(u)) z /dr versus a (left) and i?o _1 
(right) in the bulk. For co-rotation and only slight counter-rotation, there is a nearly 
linear relationship linking Ro~ l (or a) with d{uj) z /df. A black straight line has been 
added in this region to artificially extend this relationship up to d(ui) z /df = 0. There is 
also a plateau of d(ui) z /dr at counter-rotation. Here the r-slopes are smaller, i. e., the 
profiles flatter, for increasing Ta, reflecting an increased convective transport. 



We now come back to the corotating regime a < and want to connect the numerically 
found approximately linear relationship between the slope of w(r) in the bulk and Ro^ 1 
with the dynamical equation (|2.2I) . which for the (^-component of the velocity can be 
rewritten as 

dtue + u r d r ue H — -dgug H — — - + u z d z ug = Ro^ 1 u r — dgp . (6-2) 
r r 

The linear relationship can be obtained if we assume that the Coriolis force term i?o _1 u r 
and convective term u r d r ug +u r ug/r balance each other, i.e., we assume that the axial, 
azimuthal, and temporal dependences are small in eq. (|6.2I) . which then boils down to 
u r (d r ug + ug/r) ~ u r i?o _1 . Next, we use the fact the radial velocity component u r - the 
wind - in its non-dimensio nal form is cons t ant al ong the present Ta-range (cf. Section 01 
seen also in experiment of Huisman et all ( 20121) ). Therefore, and as ug hardly depends 
on i?o _1 , an increased Coriolis force can only be balanced by a larger slope d r ug. The 
only alternative is that the wind u r vanishes altogether, u r = 0, and the flow state returns 
to the purely azimuthal, laminar case. 

To further substantiate this argument, we now decompose the flow field into a {t, 9, z)- 
averaged mean azimuthal flow component Ug, depending on the radial position r only, 
plus fluctuations u' , as well as a decomposition of the pressure into a mean P plus the 
pressure fluctuations p' . Inserting these Reynolds type decompositions into the radial 
and azimuthal momentum equations - in which besides Ug only its r-derivative survives 
and ignoring viscosity for now, we arrive at the following equations: 

?/' U 2 Uaii' a' 2 
d t u' r +u'd r u' r +^-d e u' r °- e -+u'd z u' = -driP+p^-Ro-^Ug+u'g), (6.3) 

r r r r 



d t u' g + u' r d r Ug + u' r d r u' g + ^dgu'g + ^ + ^ + u'd z u' e = Ro- l u' r - dgp'. (6.4) 

r r r 

It is important to note that U r and U z are both equal to zero, so u r = u' r and u z — u' z . 
As long as i?o _1 is larger than Ro~^ t , we assume that already the mean flow contributions 
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alone balance in eqs. (|6.3|) and (|6.4|) . 

n 2 
r 

and 



-d r P + Ro^Ue , (6.5) 



u' r drUe = u r + ~ -< (W 1 + ^) . (6.6) 

As in the bulk r ~ r a is almost constant, the linear relationship 3 r CL> oc const — i?o _1 
between Ro" 1 and bulk slope d(u)) z /dr is obtained. 

Figure I21bl displaying this linear relationship, can be used to obtain a quantitative 
estimate for optimal transport for large Ta. We can see two distinct features in the slope 
d(u)) z /dr versus Ro^ 1 curve. There is a plateau, where the value of the slope is linked 
to Ta (and therefore to the viscous term in the equation of motion), and there is a line 
to the right of the plateau where the value of d(u>) z /dr is independent of Ta and thus 
linked only to the Coriolis force. From the previous discussion and from the experimental 



evidence of Ivan Gils et al\ (|2012f ) we know that optimal transport is linked to the flattest 



w-profile. We can interpret the shift of i?o~ p 1 t with Ta as that value of Ro^ 1 where the 
plateau meets the co-rotation linear relationship, i.e. the flattest possible w-profile that 
does not break the large scale balanced discuss before. If Ro^ 1 becomes more negative, 
the profile would have to become flatter to keep on satisfying the large scale balance. As 
this does not happen, the transport decreases for more negative Ro^ 1 . 

With increasing Ta, the value of d{ui) z /df at the plateau increases, and the curves cross 
at a smaller value of Ro , which corresponds to a shifted maximum. Eventually, the 



platea u value of d(oj) z /df will tend to zero as seen in the experiments of Ivan Gils et 



<|2012r ) . and the co- rotation line will cross the plateau at the x-axis. We can extend 
the straight line to get an estimate for when this happens and obtain Ro~ pt (Ta — > 
co) = —0.22, corresponding t o a ov t ~ 0.38, not far from the experimental value a op t ~ 



0.33 ± 0.05 of Ivan Gils et ali ((2012) 



If i?o _1 is too negative, Ro^ 1 < Ro~^} t , the large scale balance of equation (|6.6p can 
no longer be satisfied. This can be seen as the Coriolis force now has values which would 
require a smaller (or even a negative) value of dui /dr for the balance to hold. Since this 
is not possible to accommodate, bursts will originate from the outer cylinder towards 
the inner cylinder, because the flow tries to accommodate a large Coriolis force. These 
bursts increase in importance until they end up stabilizing parts of the flow, or even the 
whole flow which will drastically reduces the transport. Therefore, a maximum transport 
is reached just when the Coriolis force balances the large scale convective term. If it is 
further increased, stabilized regions start appearing. This is linked to the appearance of 
a neutral surface, which is analyzed in the next section. 



6.5. Neutral surface 

In this subsection, we will take a break from using the rotating reference frame, and 
return to the inertial reference frame to analyze the neutral surface which is defined as 
that surface where, in the inertial reference frame, u) is zero. This surface only exists 
for non-negative values of a and coincides with the outer cylinder in the case of pure 
inner cylinder rotation. In an inertial reference frame, it marks the division between the 
Rayleigh (inviscid) stable and unstable regions. This means that this surface separates 
two regions, an unstable inner region and a stable outer region. In the stable region, 
perturbations to the azimuthal flow (both large scale wind as in Taylor vortices and 
small scale perturbations such as plumes) cannot grow. Therefore we expect this surface 
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to play a significant role for the beh avior of the flow . It wa s already shown to be important 
in controlling optimal transport by van Gils et al. ( 20121 ). 

In general, the position of the neutral surface depends on the height z. Figure |2"21 
shows contour plots of the angular velocity with tn{z) indicated. Indeed, rjv(z) shows a 
strong axial dependence, showing heights with positive or negative u r , at which the neu- 
tral surface is pushed more outside or inside, respectively. This strong axial dependence 
of r/v is a measure for the vortex strength and becomes weaker, when the vortices lose 
importance at very high Ta. By comparing Figures [52] and [M] the effect of the Taylor 
number on the position of the neutral surface can be seen. Its distortion happens at larger 
a for increasing Ta, as expected. 

If a is large enough, the vortices are no longer able to penetrate the whole gap. There 
the neutral surface separates the Rayleigh stable and unstable regions. The vortices are 
mainly located in the unstable range, but enter partially also into the Rayleigh stable 
region. Being restricted to part of the gap, they also shrink in horizontal direction. 
Because the vortices try to remain as square-like as possible, their height (wave length) 
also shrinks, allowing new vortices to appear in the available given height. This is visible 
in the right panel of figure [231 These vortices are also associated with a stronger wind. If 
the value of a is not very high, they thus will fill up again the distance between the two 
cylinders but with a distorted aspect ratio. A zoom-in of this effect can be seen in the 
middle panel of figure 1231 This causes both the rise in Nu u for positive a seen in Figure 
[13] around a « 0.3 at low Ta and the crossovers seen in Figure [25a] 

Next, in addition to the temporal and azimuthal average, we also average tn{z) in 
axial direction and call this average r/v- Figure [25a] shows how fjv varies with a and Ta. 
The position of the neutral surface for the laminar, purely azimuthal flow is plotted for 
comparison. For slight counter-rotation and fixed a, the mean neutral surface is increas- 
ingly pushed towards the outer cylinder with increasing Ta due to enhanced turbulence. 
On the other hand, with increasing a the Coriolis force pushes the neutral surface more 
and more towards the inner cylinder. Once the neutral surface reaches the laminar and 
purely azimuthal flow value, the flow is stabilized. 

The curves fjy(a) for different Ta can also cross. At a constant rotation ratio, some of 
the lower Ta have a neutral surface which is further away from the inner cylinder than 
for some of the higher Ta. This is due to changes in the number and strengths of vortex 
pairs in the flow, which happen earlier for lower Ta. With a further increase of a the 
trend reverses again, since respectively smaller values of a stabilize the flow at already 
lower Ta. In an inertial reference frame, this simply means that there is no longer a radial 
velocity which can push the neutral surface outwards, so r/v falls back to the laminar, 
purely azimuthal flow value. 



7. Conclusions 

An extensive DNS exploration of the parameter space of a Taylor-Couette system 
at Taylor numbers in the range of 10 4 < Ta < 10 8 was presented. First the code was 
validated versus existing numerical and experimental data. After this, the transition from 
the laminar but still purely azimuthal regime to the Taylor vortex state was analyzed. 
The regime where these vortices dominate the flow was studied in detail, revealing scaling 
laws between the Taylor number Ta, the angular velocity flux Nuu, — 1, and the wind 
Reynolds number Re w . These scaling laws ceased to be valid when the Taylor number 
was increased beyond Ta 3 • 10 6 . At this driving strength the coherence structures 
become so small that they lose importance and are no longer the dominating feature of 
the flow. 
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Figure 22: Contour plots of the dimensionless angular velocity U)(r, z) with indicated 
neutral surface (black line) for Ta = 1 • 10 5 . Left: a = 0.2 and three vortex pairs. Middle: 
a = 0.4 and four vortex pairs with a non-square aspect ratio and a highly distorted 
neutral surface. Right: a = 0.6, four vortex pairs which do not penetrate the whole gap. 
See also figure [23] for a zoom-in. 




Figure 23: Zoomed-in contour plots of ui(r, z) for Ta = 10 5 with the neutral surface 
indicated as a black line. Left: a = 0.2, normal state with three vortices. Middle: a = 0.4, 
the distorted vortices are strong enough to fully penetrate the gap. Right: a — 0.6, the 
distorted vortices cannot penetrate the whole gap due to the stabilizing effects beyond 
the neutral line. Moreover, with increasing a the distance between the vortex centers 
shrinks. 
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Figure 24: Contour plots of cu(r, z) with neutral surface indicated for Ta = 4 • 10 6 and 
a-values (left to right) of 0.2, 0.4, and 0.6. All of them have three vortex pairs. 
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Figure 25: Left panel, the neutral surface location f^v of the time, azimuthally, and axially 
averaged angular velocity versus rotation ratio a for various Ta. The neutral radius vn 
is moving inwards with increasing counter- rotation (increasing a), but is pushed back to 
the towards the outer cylinder for increasing Ta at given a. The black line shows the 
position of neutral radius r/v.o in laminar, purely azimuthal flow. Right panel, fjv — r/v.o 
is shown, quantifying the pushing of the neutral line towards the outer cylinder with 
increasing Ta. 
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Then the effect of the outer cylinder rotation on these scaling laws was analyzed. If 
both cylinders are co-rotating, the scaling laws were (slightly) modified, but for counter- 
rotating cylinders no significant differences could be seen. After the shrinking of the 
coherent structure and loss of their importance the value for optimal transport a op t 
shifted towards counter-rotation. This drift is expected to continue at higher Taylor 
numbers and will be the course of future DNS investigations. 

Next, the behavior of local flow variables was studied. Analyzing the profiles Lu(r) sheds 
light on the two transport mechanisms, convective and diffusive, cf. the two contributions 
in (|1.4|) . The optimal transport of uj could be linked to a balance between the Coriolis 
force and the inertial terms in the equations of motion. This balance is best achieved 
when the bulk profile is flattest. 

The outer boundary layer of the w-profile is much thicker than the inner boundary 
layer. The quality of the approximation of the w-profiles by three straight lines was 
found to improve with increasing Ta, as the (bulk-)turbulence becomes stronger. But 
although the bulk is turbulent, the boundary layers are still of Prandtl-Blasius type. 
TC flow only reaches the ultimate state, if also the boundary layers undergo a shear 
instability and become turbulent, too. The present analysis showed that this transition 
is expected to happen in the Ta range between 10 8 and 10 9 , which is just outside the 
range of the present DNS. It will be analyzed in future work. 
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